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Abstract 

In this paper we present a uniform method for translating an arbitrary 
nondeterministic finite automaton (NFA) into a deterministic mass action bio¬ 
chemical reaction network (BRN) that simulates it. The BRN receives its input 
as a continuous time signal consisting of concentrations of chemical species that 
vary to represent the NFA’s input string in a natural way. The BRN exploits 
the inherent parallelism of chemical kinetics to simulate the NFA in real time 
with a number of chemical species that is linear in the number of states of the 
NFA. We prove that the simulation is correct and that it is robust with respect 
to perturbations of the input signal, the initial concentrations of species, the 
output (decision), and the rate constants of the reactions of the BRN. 


1 Introduction 

Molecular programming combines computer science principles with the information 
processing capabilities of DNA and other biomolecules in order to control the struc¬ 
ture and behavior of matter at the nanoscale. Molecular programming has its origin 
in Seeman’s development of DNA nanotechnology in the 1980s [28] (indeed, “molec¬ 
ular programming” and “DNA nanotechnology” are still nearly synonymous), but 
the held has made progress in the present century at a rate whose increase is rem¬ 
iniscent of Moore’s law. The achievements of molecular programming are far too 
numerous to survey here, but they include the self-assembly of virtually any two- or 
three-dimensional nanoscale structure that one wants to prescribe [STHHIIHlEnilS], 
DNA strand displacement networks that simulate logic circuits and neural networks 
pn |23[ [25]. and molecular robots that perform various functions while either walking 
on nanoscale tracks or boating free in solution [221 ESI EHl El 1131 El EZ]- All this 
has been achieved in real laboratory experiments, and applications to synthetic biol¬ 
ogy, medicine, and computer electronics are envisioned. Theoretical progress includes 
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demonstrations that various molecular programming paradigms are, in principle, Tur¬ 
ing universal [351 [H EDI Uni Ha [36], thereby indicating that the full generality and 
creativity of algorithmic computation may be deployed in molecular and biological 
arenas. 

Our objective in this paper is to begin mitigating the “in principle” of the pre¬ 
ceding sentence. This is important for two reasons. First, although such theoretical 
results are steps in the right direction, processes that require unrealistically precise 
control of unrealistically large numbers of molecules simply cannot be implemented. 
Second, processes that can be implemented, but only with inordinately precise con¬ 
trol of parameters are inherently unreliable and hence inherently unsafe in many 
envisioned applications. Our objective here is thus to identify a class of computa¬ 
tions that can be implemented robustly in the molecular world, i.e., implemented in 
such a way that they will provably perform correctly, even when crucial parameters 
are perturbed by small amounts. Future research can then strive to enhance this 
robustness and to extend the class of computations that enjoy it. 

In this paper we give a uniform method for translating nondeterministic hnite 
automata to biochemical reaction networks that implement them robustly. Nonde¬ 
terministic hnite automata (NFAs) are over half a century old [26| and far from Tur¬ 
ing universal, but they have many applications and remain an active research topic 
Hiiniis]. Applications of NFAs that are likely to extend to molecular programming 
include their uses in monitoring and parsing large data streams and in implementing 
and verifying secure and/or safe communication protocols. Biochemical reaction net¬ 
works (BRNs) are also at least half a century old [3|. Originally-and still often-called 
chemical reaction networks (CRNs), their role in molecular programming did not be¬ 
come fully apparent until recently, when Soloveichik, Seelig, and Winfree [31] showed 
that there is a systematic method for translating an arbitrary BRN, which is an ab¬ 
stract mathematical object, into a set of DNA strands and complexes that simulates 
the BRN via toehold-mediated strand displacement. This method has been rehned 
and formulated as a compiler [8], and BRNs are now the programming language of 
choice for many molecular programming investigations. (The fact that this systematic 
implement at ion-as opposed to earlier ad hoc implementations of individual reaction 
networks-uses DNA has probably contributed to the CRN-to-BRN terminology shift 
that seems to be underway.) The two most widely used semantics (operational mean¬ 
ings) for BRNs are deterministic mass action semantics and stochastic mass action 
semantics. In this paper we use deterministic mass action, which implies that the 
state of a BRN at any time is determined by the real-valued concentrations of its 
molecular species at that time. 

An NFA is a real-time device that reads its input string sequentially, left to right, 
changing states appropriately in response to each symbol prior to reading the next 
symbol. Accordingly, we translate each NFA to a BRN that receives the NFA’s input 
string formatted as a continuous time concentration signal consisting of concentrations 
of input species that vary to represent the input string in a natural way. Using the 
inherent parallelism of chemical kinetics, our BRN implements the NFA in real time, 
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processing each input symbol before the next one arrives, and it does so with a 
number of molecular species that is linear in the number of states of the NFA that 
it implements. Specihcally, if the NFA has q states, s symbols, and d transitions, 
then our BRN has 4g + s + 2 species and 5g + d reactions. The compiler of [8] 
would then translate this into a DNA strand displacement network consisting of 
28g + 4d + 2s + 6 single strands of DNA. Our translation thus appears to make small 
NFA implementable in laboratories now and NFAs of modest size implementable in 
the near future. 

Most importantly, our BRN’s correct implementation of the NFA is robust with 
respect to small perturbations of four things, namely, its input signal, the initial 
concentrations of its species, its output signal, and the rate constants of its reactions 
(acceptance or rejection of its input string). One key to achieving this robustness is 
a memory-refresh technique akin to the approximate majority algorithm of PdE]. 

The NFA and BRN models and their semantics are reviewed in detail in section 
2. The robustness properties of our BRNs are explained in section 3. Our translation 
from NFAs to BRNs is specihed in section 4, and section 5 gives the proof that our 
BRNs are robustly correct. 


2 Review of NFAs and BRNs 

Finite automata are ubiquitus in computer science, but details and notation vary 
from context to context, so we briefly review the specihc model used in this paper. 
(See, e.g., |2l].) 

A nondeterministic finite automaton (NFA) is an ordered 5-tuple 

M = (Q,S,A,J,F), 

where Q is a hnite set of states] S is a hnite input alphabet] / C Q is the set of initial 
states] F C Q is the set of accepting states] and A : Q x S —)■ V{Q) is the transition 
function. Here we are using the notation V{Q) for the power set of Q, i.e., the set 
of all subsets of Q. When convenient we identify the transition function A with the 
ternary relation 

A = {(g,a,r) eQxExQ|rG A(q',a)}, 

which is the set of all transitions of M. Informally, the size of M is determined by 
the three cardinalities |Q|, |S|, and |A|. 

The extended transition function of the above NFA M is the function A : V{Q) x 
S* —)■ V{Q) dehned by the recursion 

A(A,A) = A, 

A{A,wa) = A(g, a) 

qGA{A,w) 
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for all A C Q, w G S*, and a G S, where A is the empty string. The NFA M accepts 
an input string ta G S* if A(J, ta) flF 7 ^ 0, i.e., if there is a chain of transitions leading 
from some state in I to some state in F. 

Biochemical reaction networks are not as widely known in computer science, so 
we review them less tersely. 

Formally, we £x a countably infinite set S = {Xq, Xi, X 2 ,...}, whose elements 
we call species. Informally, we regard each species as an abstract name of a type of 
molecule, and we avoid excessive subscripts by writing elements of S in other ways, 
e.g., X, X, Z, X, X, etc. 

A reaction over a finite set S' C S is formally a triple 

p = (r, p, /c) G X X (0, cx)), 

where is the set of functions from S into N = {0, 1 , 2 ,...}, and r 7 ^ p. Since S is 
hnite, it is natural to also regard elements of as vectors. Given such a reaction p, 
we write r(p) = r, p(p) = p, and k{p) = k, and we call these three things the reactant 
vector, the product vector, and the rate constant, respectively, of the reaction p. The 
species in the support set supp(r) = {X G S'|r(X) > 0} are the reactants of p, and 
the species in supp(p) are the products of p. 

We usually write reactions in a more intuitive, chemistry-like notation. For ex¬ 
ample, if S' = {X, Y, Z}, then we write 

X + Z ^2Y + Z 

for the reaction (r, p, k), where r, p : S' —)■ N are defined by 

r(X) = r(Z) = l, r(X) = 0, 
p(X)=0, p(X) = 2, p{Z) = l. 

The net effect of a reaction p is the (nonzero) vector Ap = p(p) — r(p) G 
A species Z satisfying r(p)(Z) = p(p)(Z) > 0, as in the example above, is called a 
catalyst of the reaction p. 

A biochemical reaction network (BRN) is an ordered pair 

N = {S, R) 

where S' C S is finite, and i? is a finite set of reactions over S'. 

There are various semantics (i.e., operational meanings) that can be usefully as¬ 
cribed to a BRN N = (S', R). In the deterministic mass action semantics (also called 
deterministic mass action kinetics) that we use in this paper, a state of the BRN N 
is a vector x G [0,cx))‘^. For each X G S' we call the nonnegative real number x(X) 
the concentration of the species X in the state x. 

In deterministic mass action semantics, the state of a BRN, and hence the con¬ 
centrations of its species, vary continuously over time. This makes the following 
notational convention convenient. Concentrations of species X, X, etc., are denoted 
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by the corresponding lower-case letters x, y, etc. Thus, for example, we write y{t) 
for the concentration of Y at time t. This convention allows us to write a state of a 
BRN N = (S', R) as a vector x = (y | F G S) G [0, ooY. 

Let N = {S,R) be a BRN. For each p E R and x G [0,oo)'^, the {deterministic 
mass action) rate of p in state x is 

ratex(p) = k{p)x''^f’\ (2.1) 


where we write 

Xr(p) = 

Yes 


k 

For example, if p is the reaction X R Z ^ 2Y + Z, then ratex(p) = /i;x(X)x(Z). 
Intuitively, the frequency with which an X and a Z encounter one another is propor¬ 
tional to x(X)x(Z), and the constant of proportionality k summarizes other factors, 
not depending on x (e.g., temperature, salinity of solution, properties of X and F), 
that also govern the rate at which p occurs. 

For each species Y E S dehne the {deterministic) mass action function Fy ■ 
[0, cx))‘^ -)■ M by 

Fy(x) = y^ratex(p)Ap(y) (2.2) 

peR 


for all X G [0, oo)'^. Then FV(x) is the total rate at which the concentration of Y 
is changing in state x. Hence, if at time t the BRN N is in state x(t), then the 
concentration y{t) of each species Y must obey the ordinary differential equation 
(ODE) 

y'{t)=FY{^{t)). (2.3) 


If we let Sy be the ODE (2.3) for each Y E S, then the {deterministic) mass action 
system of N is the coupled system 


{Sy\Y E S) 

of ODEs. If we dehne the vector-valued function F : [0, cxd)*^ —)■ by 

F(x) = (Fy(x) \YeS) 


(2.4) 


( 2 . 6 ) 


for all X G [0, cx))*^, then the mass action system (2.4) can also be written in the vector 
form 

x'(t) = F(x(t)). (2.6) 


It is clear by inspection of (2.1)-(2.6) that the mass action system (2.6) is 


autonomous, meaning that its right-hand side only depends upon the time t indirectly, 
via the concentrations y{t), and (ii) polynomial, meaning that its right-hand side is a 
vector of functions that are polynomial in their arguments y{t). 

In deterministic mass action semantics a BRN N is initialized to some state at time 


0, and this state then evolves according to the mass action system (2.6). Accordingly, 
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we define an initialized biochemical reaction network (IBRN) to be an ordered pair 
(iV, xo), where iV is a BRN and xq is a state of N. The {deterministic) mass action 
initial value problem {mass action IVP) of an IBRN (iV, xq) is the initial value problem 
(IVP) consisting of the mass action system (2.6) of N together with the initial value 
condition 

j/(0) = xo(R) for each Y e S. (2.7) 


By the standard existence-uniqueness theory for ODEs EES, this mass action IVP 
has a solution x(t) that is defined for all t G [0,6) for some b G (0, cxd], and this 
solution is unique. It is not difficult to show, then, that x(t) G [0, cxd)'^ holds for all 
t G [0,6), i.e., that concentrations remain nonnegative. In deterministic mass action 
semantics, this solution x : [0, 6) —)■ [0, oo)^ is the operational meaning of the IBRN 
{N, Xq), and the family of such solutions, as xq ranges over [0, oo)^, is the operational 
meaning of the BRN N. The BRNs defined in this paper are all very well behaved, 
so that 6 = oo, i.e., x(t) is well defined for all t G [0, oo) and all initial values xq 
considered here. 

Further discussions of biochemical reaction networks with deterministic mass ac¬ 
tion semantics appear in [121 HSl El E] • 


3 Robustness 

The objective of this paper is to translate an arbitrary NFA M = (Q, S,A,J, F) 
into a BRN N = (S', R) that is robust with respect to four types of perturbations, 
namely, perturbations of the input signal, the initial concentrations of species, the 
output (decision) of the BRN, and the rate constants of reactions. The present section 
describes these types of robustness separately. 

3.1 The input signal. An NFA M = {Q,T,, A, I, F) receives its input 
w G S* from an external agency (user) that is not part of the NFA model. This 
input w is presented to M sequentially, reading left-to-right, with M changing states 
appropriately in response to each symbol of w before reading the next symbol. Ac¬ 
cordingly, the BRN N = {S, R) that simulates M receives the input string tc G S* 
formatted by an external agency as a continuous time signal consisting of concentra¬ 
tions of input species that vary to represent tc in a natural way. This input signal and 
the perturbations of it with respect to which our BRN is robust require some detailed 
explanation. 

It is the user’s responsibility to provide the input signal, i.e., to control the con¬ 
centrations of the input species. This, in turn, implies that our BRN should not 
interfere with the user’s control. Hence the input species only appear as catalysts of 
the reactions in our BRN. 

Our BRN N has |S| -|-2 input species, namely, a species Xa for each input symbol 
a G S and two special input species, Xr (“reset”) and Xc (“copy”). Intuitively and 
idealistically, N would receive an input string w = UqUi ■ ■ ■ a„_i, where each G S, as 
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a 0/1-valued “square-wave” (i.e., step-function) representation of the “padded” string 
X = aocraicr ■ ■ ■ an-icr, where the square-wave representation of x is as depicted in 
Figure 


x{t) 




Figure 1: Idealized square-wave format of padded version w = raocraic- ■ -ran-ic of 
input string w = UqUi ■ a„_i G S*. All input species concentrations are 0 at all times 
when the hgure does not indicate that they are 1 . 

There are two things wrong with this idealization. First, step functions are not 
continuous, and it is unrealistic to hypothesize an agency that varies concentrations 
discontinuously. More importantly, it is unrealistic to hypothesize an external agency 
that varies concentrations according to any scheme that does not tolerate small errors. 
We thus allow the following relaxed input format. An input signal is an assignment 
of a continuous concentration function x : [ 0 , oo) —?• [ 0 ,oo) to each of the |S| -|- 2 
input species X. In the following we £x a phase duration r > 0. Let e G (0, ^). An 
input signal is an e-admissible encoding of an input string w = aoUi ■ ■ ■ Un-i for M if 
it satisfies the following conditions for all k gN. (See Figure [^) 

(1) All input species have concentrations less than 1 -|- e at all times. 

(2) All input species have concentrations less than e at time kr. 

(3) All but at most one input species X have the property that, at all times t G 
[kr, {k + l)r], the concentration x{t) is less than e. Such species are not present 
during phase k. 

(4) If any input species is an exception to (3), then its concentration exceeds 1 — e 
throughout the time interval [kr + ^, kr + ^]. 

(5) If fc = 3i, where 0 < i < n, then input species Xr is present during phase k. 

( 6 ) If A; = 3i + 1, where 0 < i < n, then input species is present during phase 

k. 

(7) If k = 3i + 2, where 0 < i < n, then input species Xc is present during phase k. 

( 8 ) If A; > 3n, then no input species is present during phase k. 
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kr [k + l)r kr [k + l)r 


(a) 


(b) 


Figure 2: e-admissible input signals. In (a) the input species X is present during 
phase k. In (b) the input species X is not present during phase k. In either case, any 
concentration signal conhned to the interior of the white region is e-admissible 


3.2 Initial concentrations. Given two initial states xo,Xq g [0,cxd)'^ of a 
BRN N = (S', R), we dehne the distance between these to be 

\\^*o - xoll = max |x*(F) - xo(F)|. 

Intuitively, this is the “amount by which xq has been perturbed” if xq is replaced by 


3.3 Output. As seen in subsection 3.1 above, our BRN N receives an input 
X = a^ai ■ ■ ■ Qn-i G F* for M as an e-admissible signal that arrives in Sn phases, 
each of duration r. After this time, no input species is present, i.e., all input species 
have concentrations less than e. At any time after one additional phase, i.e., any 
t > (3n l)r, a user of our BRN N may observe one or more of the concentrations 
yqit) of those species Yg that correspond to states q E Q of the NFA M. This 
observation itself is subject to error, i.e., it may be a value yq{t) that differs from 
yq(t). Let 7] G (0, b). We say that a state observation scheme yq h->■ jjq is rj-admissihle 
if, for every t > (3n -|- l)r and every state q E Q of M, the observed concentration 
yq(t) satishes 


IVgit) - yqit)\ < V- 














3.4 Rate constants. In order to state the sense in which our BRN’s simula¬ 
tion of an NFA is robust with respect to perturbations of rate constants, we need a 
slight generalization of the BRN definition. 

A time-dependent reaction over a hnite set S' C S is formally a triple 

p = (r, p, k*) G X X C([0, cxd), (0, cxd)), 

where (^([O, oo), (0, oo)) is the set of continuous functions from [0, oo) to (0, oo). As 
before, we write r(p) = r, p(p) = p, and k*{p) = k*, and we use more intuitive 
notions like 

X + Z ^2Y + Z, 

remembering that k* is now a function of time, rather than a constant. A time- 
dependent BRN (tdBRN) is then an ordered pair 

N* = (S,R*), 


where S C S is finite, and R* is a finite set of time-dependent reactions over S. The 
deterministic mass action semantics of a tdBRN N* = (S, R*) are then defined in the 
obvious way, rewriting (2.1)- (2^ 


as 


ratex(p)(t) = k*{p){t)xl^^f’\ 

(3.1) 

ratex(p)(t)Ap(y'), 

p£R* 

(3.2) 

y\t) = Fy(x(f),f), 

(3.3) 

(Ty 1 A G 5), 

(3.4) 

F(x,t) = (Fy(x,f) |y G^), 

(3.5) 

x'{t) = F(x(f),f). 

(3.6) 


Note that (3.6) is not autonomous. 

Let N = (S, R) be a BRN, and let 5 > 0 be a positive real number. A 6- 
perturbation of A^ is a tdBRN N* = [S, R*) in which R* is exactly like R, except that 
each reaction (r, p, k) is replaced by a time-dependent reaction (r, p, k*) satisfying 


\k*{t) -k\<6 


(3.7) 


in which an adversary is allowed 


for all t G [0, cx)). 

Intuitively, a (5-perturbation of A^ is a tdBRN N* 
to vary each rate constant k, subject to the constraint (3.7). 

Putting this all together, let e, 77 G (0, \) and (5 > 0; let M = (Q, S, A, /, F) be an 
NFA; and let (A, xq) be an IBRN, where N = {S,R). We say that (A, xq) {e,p,6)- 
robustly simulates M if the following conditions hold. (Recall that we have fixed a 
phase duration r > 0.) 

(i) A has the |E| -I- 2 input species described in subsection 3.1, and these only 
appear as catalysts of reactions in R. 
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(ii) There is a state species E S for each state q E Q. 

(iii) For each q E I, ^o(Yq) = 1. For each q E Q\I, xo(Fg) = 0. 

(iv) For every (5-perturbation N* of N, every input string w E T,* for M, every e- 
admissible encoding of w as an input signal, every initial state Xq of N* satisfying 
||xq — xqII < e, every T^-admissible state observation scheme i/q h-)■ ijq, and every 
time t > (3n -|- l)r, the observed concentrations jjqit) in (xq, N*) satisfy 

2 

Vqi^) > 3 e 

and 

yq{t) < 3 <1 i 

4 Translating NFAs to BRNs 

We now describe our method for translating an arbitrary NFA M = (Q, S,A,/, F) 
into a BRN N = {S,R). 

4.1 Species Our translation uses exactly 4IQI-I-|S|-|-2 species. Of these species, 
there are four species types that are useful to distinguish. 

1. For each state, q E Q, there is a state species Yq whose presence in the BRN 
corresponds to the NFA being in state q. 

2. For each state, q E Q, there is a portal species, Zq, which is used as a buffer to 
facilitate state transitions in the BRN. 

3. In order to support e-admissible encodings of input strings, there are exactly 
|S| -|- 2 input species. There is a species Xa for each symbol a G S, and the two 
special input species Xr and Xf.. The species Xr intuitively is used to “reset” 
the BRN in order to prepare to process the next input symbol. The species X^ 
is used along with the portal species in the BRN to “copy” the current value of 
each portal species, Zq, into the corresponding state species, Yq. 

4. For each state species Yq and portal species Zq, there are dual species Yq and Zq. 
We refer to the species Yq, Zq as basic species in order to further distiguish them 
from their duals Yq, Zq. Intuitively, a dual of a basic species is one that has 
exactly the opposite operational meaning, i.e., when Yq has high concentration, 
Yq has low concentration and vice versa. 
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4.2 Reactions Our translation uses exactly 5|Q| + |A| reactions and they are 
enumerated below. 


• Vg G Q 

Xr + Zg%Xr + Zg (4.1) 

• V(s, a,q) e A 

X^ + Y, + Z^^Xa + Ys + Zg (4.2) 

• 'iq & Q 

X, + Zq + Yq^X, + Z, + Yg (4.3) 

X, + Zg + Yg%X,+ Zg + Yg (4.4) 


'^q e Q 

2Yg + Yg 

2Y, + K 


A14 


^4 


31', 

31 ; 


(4.5) 

(4.6) 


Ideally, our BRN would perform three different computations in sequence. 


1. Reset the portal species to 0 (using reactions (4.1)) 

2 . Compute the state transition function and store the result in the portal species 
concentrations (using reactions (4.2)) 


3. Copy the value of the portal species into the state species (using reactions (4.3)- 

0) 

In reality, all the reactions are active simultaneously and cause errors to accumulate 
to catastrophic levels. Hence the reactions (4.5)-(4.6) are included to mitigate these 
errors. These reactions compute an approximate majority algorithm to keep the 
concentrations of the state species and their duals at appropriate distances. 


4.3 Example. Finally we give a small example to demonstrate the transfor¬ 
mation from a small NFA that decides if a binary string’s second to last bit is a 
1. (See Figure i) The NFA has 3 states and 5 transitions, so the BRN contains 
4-4-1-2-1-2 = 16 species, and 5 ■ 3-|-5 = 20 reactions, and is shown in Figure]^ below. 
Note the species are not listed explicitly but are implicitly listed within the reactions. 
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Figure 3: Example NFA with 3 states and 5 transitions that decides if the second 
to the last bit is a 1. 



Figure 4: BRN reactions for example NFA. 
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5 Robust correctness of translation 


We now give a complete proof that for an arbitrary NFA M our BRN translation 
N robustly simulates the NFA, provided some modest constraints are placed on the 
parameters. We hrst state the main theorem here with proof using lemmas that 
follow. 

Theorem 5.1. Let e ,?7 G (0, |), let r, 5 > 0, and let M = (Q,S, A,/, F) be an NFA. 
Let N = [S, R) be the BRN constructed from M in section and let xq be a state 
with 


^qeQ\I, 

Vg G Q, 


Xo(>"g) = 1, 
Xo(Rg) = 0, 
Xo(Zg) = 0, 


xo(^g) = 0, 
Xo(L;) = 1, 
Xo(Zg) = 1. 


If the constraints from lemmas 5.27| and |5.28| 
bounded such that 1 — 7 > | + 7 and 7 < | 
simulates M. 


hold, and the helper constant 7 is 
— 1 ] then (A, xq) (e, 7 , (5)-robustly 


Proof. Assume the hypothesis. By lemmas 5.27 and 5. 


we have that for all tc G S* 

the inductive hypothesis (dehned below) $(ic) = true. Since 1 — 7 > | + 7 and 


7 < 


1 


rj, the theorem holds. 


□ 


The rest of this section is devoted to proving the lemmas required for theorem 5.1 
Let N* be a (5-perturbation of N, and let Xg be a state such that ||xq — xo|| < e. 
Let w G S* be an input string of M, and let the input species of N conform to an 
e-admissible encoding of w. 

Recall that the tdIBRN (A*,Xg) specihes an IVP that has a unique solution. It 
is easy to verify that the set of ODEs generated by N* are as follows. 


dyq 

dt 

dZq 

dt 

dZq 

dt 

dyq 

dt 


klXc [Zqyq - Zqyq] + [|/g - ^/g] 


klXaysZq 

{s,a,q)eA 


dZq 

dt 

dyq 

dt 


k^XrZ 


q 


(5.1) 

(5.2) 

(5.3) 

(5.4) 


In order to prove that N (e, rj, (5)-robustly simulates M, it suffices to show that 
for all t > (3|w| -|- l)r, the solution to (A*,Xg) satishes 


1 5 2 

VqW > 3 + ^’ 

yqit) < ^ 


Vg G A{I,w), 
Vg ^ A{I,w), 
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(5.5) 

(5.6) 










We prove this in two steps. First, we prove by induction on tc G E* that at 
time t = 3|ta|r, the solution to (A^*,Xq) “almost” satishes the statements (5.5)-(5.6). 
Second, we prove that the approximate majority structure within N* achieves the 
necessary bounds by time t = (3|ta| + l)r and maintains them indehnitely. Together, 
these two parts trivially satisfy the statements (5.5)-(5.6). 

For clarity, we state our induction hypothesis as the Boolean function $ : E* —>■ 
{true, false}. Let 7 G (0,^) be a constant. Dehne <h(w) = true if and only if the 
following two conditions hold. 


VgGA(/,'u;), ?/g(3|w;|r) > 1 - 7 , (5.7) 

yq^A{I,w), 2 /g( 3 |w|r) < 7 . (5.8) 

Lemma 5.2. If 7 > e, *h(A) is true. 

Proof. Recall that A(/, A) = /. Therefore, it suffices to show 


Vg G /, yq{Q) >1-7, 

Vg G Q \ /, yq{ 0 ) < 7 . 

Since ||xg — xo|| and 7 > e we have 

Vg G /, yq{Q) >l-e>l -7 

and 

Vg G Q \ /, yq{0) < e < 7 

□ 

To prove that <F(tc) holds for every string tc G E*, it remains to be shown that 
for every string tc G E* and for every symbol a G E, the following implication holds. 


<F(tc) 


<F(tca) 


(5.9) 


The majority of this section is devoted to the proof of 
series of lemmas in support of this proof. 

For each basic species Yq, Zq, let c(Yq), c{Zq) be the 


(5.9). We 


begin by proving a 


following constants. 


c(Y,) = y,(0) + j/,(0) 
c(Z,) = z,( 0 ) + z,( 0 ) 


(5.10) 

(5.11) 


Lemma 5.3. For all t G [0, cx)) and for all g G Q 


(a) yq{t)+yq{t) = c{Yq) 

(b) Zqit) + Zqit) = c(Zq) 
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Proof. Recall the ODEs for yq,yq,Zq and Zq have the following relationships. 

dyq _ _ dzq 

dt dt dt dt 

Integration and solving for the constants trivially conhrms (a) and (b) are trne. □ 
Lemma 5.4. For all q ^ Q, c{Yq), c{Zq) G [1 — e, 1 + 2e]. 

Proof. Recall every basic species Yq, Zq in the state xq is set to 0 or 1 and its Yq, Zq 
dnal is set to the appropriate 0 or 1 valne. Since ||xq — xo|| < e, it is clear that 
c(Yq), c{Zq) < 1 + 2e. Since no species can be negative, it is clear that c(Fg), c{Zq) > 
1-e. □ 


Construction 5.5. Let a > 0, 6 > 0, c > 0, and p > 0 be constants and let u(t) and 
u(t) be fnnctions dehned by 


— = uuiau — bu) — cu, 
dt ^ ^ ’ 


( 6 . 12 ) 


where u{t) = p — u{t) for all t G [0, cx)). 

2 2 

if c < 4 ^^^, then the fnnction u{t) has the following 


Lemma 5.6. In constrnction 
eqnilibrinm points. 


5.5 


^1 = 0 , 

p{a + 2b) — A 
2(a + b) ’ 
p{a + 2b) + A 


E2 = 

E3 = 


where 


2(a + b) 

A = \/p‘^a‘^ — 4c(a + b). 


(5.13) 

(5.14) 

(5.15) 


Of these three eqnilibrinm points, (5.13) and (5.15) are exponentially stable, and 
(5.14) is nnstable. 

Proof. The eqnilibrinm points of u(t) can be fonnd by solving the eqnation. 

uu[au — bu) — cn = 0 (5.16) 

It is rontine to show that the roots of ( |5.16| ) are indeed the points (5.13)-(5.15). 

Note that u{t) has three distinct real roots if and only if A is a non-zero real. The 
constant A is real and non-zero if and only if 


p^a^ > 4c(a -|- 6), 


whence. 


c < 


p^a^ 

4(a -I- b)' 
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What remains to be shown is the stability of the three points. It is well known in 
dynamical systems that if each eigenvalue of the Jacobian matrix at an equilibrium 
point has negative real part, then the point is exponentially stable. Likewise, if 
each eigenvalue has a positive real part, then it is unstable. Since the system in 
construction (5.5) is 1-dimensional, there is only one eigenvalue. 

Equilibrium point (5.13) has the eigenvalue Ai = —(c -f bp^) which is always 
negative. Therefore it is exponentially stable. 

Equilibrium point (5.14) has the eigenvalue A 2 = 2 {a+b) [p(q + 25) — Since c > 0, 
we know 0 < A < pa, therefore A 2 > which is always positive. Therefore 

it is unstable. 

Equilibrium point (5.15) has the eigenvalue A 3 = — A 2 . Since A 2 > 0, we know that 
A 3 < 0. Therefore it is exponentially stable. □ 


Corollary 5.7. For construction 5.5 and to > 0, if u{to) > E 2 , then for all f > 0, 
u{to + t) > u{to). 


Lemma 5.8. Assuming construction |5.5[ let ui,U 2 be constants such that > U 2 > 
ui > E 2 , let to > 0, and let t > 0 such that u(to) = ui and u(to + t)= U2- Then 


^ ^ TW log 


U2{Eo - Ml) 


A\Eo \_ui{Eo-U2) 




Ui{u2 - E 2 ) 
_U2{ui - E 2 ) 


(5.17) 


Proof. Recall that the ODE for u{t) is 

du 


dt 


= uuiau 


bu) 


cu 


= u{p — u) [au — b{p — m)] — cu, 


and the roots of the ODE are Ei, E 2 , and E 3 . It can be shown that 

nil 

— = (a b)u{u - E 2 ){Eo - u). 

By separation of variables and integration, we have 


J_ f _L_ 

a + bj u{u-E 2 ){Eo 



to+i 

I dt 

^0 
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Therefore, 


t = 


1 


a + 6 yi?2-^3(-^2 ~ -^ 3 ) 
1 / 1 


log 


Ez / 77' \ -^2 

U \ { u — 


u- E 2 


u 


U 2 


Ui 


a + 6 \ E 2 E^[E 2 — T/ 3 ) 


log 


M2(Mi - E 2 )\^^ (Ui{u 2 - E 3 ) 


Ui{u2-E2)J \U2{Ui-E3) 


E 2 


a + b \E 2 — E-] 


a + b 


a + b 


A 




^log 


U2(ui - E2) 
Ui{u2 - E 2 ) 

U2{Ux - E 2 ) 




^log 


E 2 \_Ui{u2-E2)\ E 3 \_U2{E3-Ui) 


Ui{u 2 - E 3 ) 

U 2 {Ui - E 3 ) 

Ui{E3 - U 2 ) 




Ui{u 2 - E 2 ) 
U 2 {Ui - E 2 ) 




U2{E3 - ^i) 

Ui{E 3 - U 2 ) 


□ 


Construction 5.9. Let a > 0, 6 > 0, c > 0, and p > 0 be constants and let u{t) and 
u{t) be fnnctions defined by 


du , 

— = uu{au — bu) + cn, 


(6.18) 


where u{t) = p — u{t) for all t G [0, cx)). 


Lemma 5.10. In constrnction 
following eqnilibrinm points. 


5.9 


2l2 

if c < then the fnnction u{t) has the 

bp - A* 


where 


_ 

^ “ 2(a + 6)’ 

_ bp + A* 

" 2(a + 6)’ 

e; = p, 


A* = y/p2^2 _ 


(5.19) 

(5.20) 

(5.21) 


Of these three eqnilibrinm points, (5.19) and (5.21) are exponentially stable, and 


(5.20) is nnstable. 


Proof. Using a similar argnment as the proof of lemma 5.6, it is easy to show that 
this lemma also holds. □ 


Corollary 5.11. For constrnction 5.9 and to ^ 0, if M(fo) < E^, then for all f > 0, 
u{to + t) < u{to). 
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Lemma 5.12. Assuming construction 5.9, let Ui,U 2 be constants such that > 
ui> U2> El, let to > 0, and let t > 0 such that u{to) = ui and u{to + t) = U2- Then 




A* \El 


U2{Ui - El) 

Ui{u2 - El) 


+ — log 
E* ^ 


Ui{El - U2) 
U2{EI - Ui) 


(5.22) 


Proof. Using a similar argument as the proof of lemma |5.8[ it is easy to show that 
this lemma also holds. □ 

Construction 5.13. Let a > 0, 5 > 0, and p > 0 be constants, and let u{t) and u{t) 
be functions dehned by 

= au — bu, (5.23) 

dt ' ^ ^ 

where u(t) = p — u(t) for all t G [0, cx)). 

Lemma 5.14. In construction |5.13| if to,t > 0, then 


u{to + t)= (1 - . 


a + b 


(5.24) 


Proof. The ODE (5.23) is trivially solvable via separation of variables and integration. 
It is routine to show that its solution matches that of equation (5.24). □ 


Corollary 5.15. In construction 5.13 
(a) if a = 0, then 


(b) if 6 = 0, then 


u{to + t) = u{to)e 
u{to + t) = M(to)e““* + p(l - e““*). 


Corollary 5.16. In construction 5.13 for to,t > 0, regardless of the initial condition 
u(to), the function u(t) is bounded in the following way. 


u{to + t) <p(^^ + e “j 


u{to + t) >p -e 

a 


—at 


(5.25) 

(5.26) 


Proof. To show (5.25), we have 


u{to + t)= u(fo)e"(“+^)* + (1 - 


a + b 


< ( 1 ) 


a + b 


< p ( e-“ + ^ 
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To show (5.26), we have 


u(to +1 ) 


= (1 - 



□ 


We now begin proving lemmas showing bounds of the solution to {N*,x.q) while 
processing input. Let d = |A|. 

Construction 5.17. Let X be an input species in {N*,x.l) that is present during 
the interval [kr, {k + l)^]. 


(6.27) 


Lemma 5.18. In construction 5.17 if X = then for all q E Q, 

Zg{{k + l)r) < 2 (^ + + ar 

where 

a = 2 d(ki + 5)e 
h={k^-d){l-e) 

Proof. Recall the ODE for Zq is 

^ ^ klXaVsZq - klXrZq. 

{s,a,q)&A 

Since d = |A|, there can be at most d reactions contributing to the sum, and by 
lemma 5.4 we know that the species Yg has at most a concentration of l+2e. Therefore, 
we have 

dz 

< dkfXail + 2e)Zq - klXrZq 

< 2dklXaZq — k^XrZq 

The constraints on the 5-perturbated rate constants ensure that for all t G [0, cx)), 
\k*{t) — k\ < 6 , whence 


< 2d{ki + 6 )XaZq 


{ks - 6)XrZq. 
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Since is present, we know that during the interval [kr + |, fcr + y] Xr > 1 — e and 
all other input species are less than e. Therefore, we have 

d z 

< 2d{ki + 5)ezq - {k^ - ^)(1 - e)zq 

= aZq — hZq 


which by corollary 5.16 yields 


^q{kT + -^) < (1 + 2e) 

Finally, during the time interval [kr + y, (fc + l)r] the ODE for Zq is bounded by 

dt - 

< a(l + 2e) 

< 2 a 


and therefore 

Zq{{k + l)r) < Zq{kT + y) + 2a 

< 2 (^ + + ar 

□ 

Construction 5.19. For q E Q, let a, b, c, and p be constants such that 

a = k^ — S, c = 2e{k2 + 5), 

h = ki + 5, p = c{Yq). 


With respect to the above constants, let Ei, E 2 , E^, and A be dehned as in lemma 5.6 


Lemma 5.20. Given construction 5.17 with X ^ Xc and construction 5.19 for q E Q, 
let yi, y 2 be constants such that E 3 > y 2 > yi > E 2 and yqikr) > yi. Then if 


then 




' 1/2 (^3 
.1/1 (^3 


2 / 2 ) 


+ 



y 1 ( 1/2 - E 2 ) 
2/2(?/i - E 2 ) 




yq{{k + l)T) > 1 / 2 - 
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Proof. Recall that the ODE for Yq is 


^ =klXc [ZqVq - ZqVq] + klVqVq [Vq - Vq] 

= klx^ZqVq - klx^ZqVq + kly^Vq - klyqf^ 

>-{k 2 + (5)e(l + 2e)yq + (^4 - - {h + 5)yqyl 

> -2e{k2 + S)yq + {k^ - 5)ylyq - {k^ + 5)yqfq 

~ VqVq V^Vq ^l/g] ^Uq- 


Let r* > 0 be a constant such that yqikr + t*) 
that 


1 / 2 . Then by lemma 5.8, we know 


T* < 


1 

A 


< T 


f \ y2{E3 

Ua [yiiEs 


yA 

2/2) 



yi{y 2 - E2) j\ 
1/2(2/1 - E2)\) 


Therefore, it is clear that yq{{k + l)r) > 2 / 2 . □ 

Construction 5.21. For q ^ Q, let a, b, c, and p be constants such that 


CL — k^ T c — 2 (/c2 T 5)6, 

b = k 4 -S, p = c{Yq). 


With respect to the above constants, let E^, E 2 , E^, and A* be dehned as in 
lemma 15.101 


Lemma 5.22. Given construction 5.17 with X ^ Xc and construction 5.21 for q G Q, 
let 2/1,2/2 be constants such that > 2/1 > 2/2 > E^ and yqikr) < 2/1. Then if 


r > 


1 

A* 



' 2/2(2/! -^1) ' 

_yiiy2 -Ei)_ 


+ 



' yijE; 


1 / 2 ) ' 

yi). 


then 

yqiik + l)T) < 2 / 2 . 

Proof. Recall that the ODE for Yq is 

^ = klXc [Zqyq - ZqPq] + klPqPq [pq - Pq] 

= k^XcZqPq - k^XcZqPq + kly^Pq - klPqfq 

< (/C 2 + 5)e(l + 2e)yq + {k^ + 5)ylyq - {k^ - S)yqfq 

< 2ik2 + 5)eyq + [k^ + 5)ylyq - [k^ - 5)yqfq 
yqVq ^//g] T ^yq- 
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Let r* > 0 be a constant such that yq{kT + r*) = y 2 . Then by lemma 5.12, we know 
that 


y 2 {yi - El) 
A* lyiiy 2 - El) 

< T 


f 4: log 


+ Tyrlog 


2/1 (^2 - 2/2) 


E*2 ^ L2/2(^2*-2/i)J 


Therefore, it is clear that yq{{k + l)r) < y 2 . 


□ 


Lemma 5.23. Given construction 5.17 with X = Xa, let 2 / 0 , and (3 be constants 
such that 


2/0 > E 2 (from construction 5.19), 
a = {ki- (5)(1 - e)2/o, 

/3 = {k^ + 5)e. 


If {s,a,q) G A and ysikr) > 2 / 0 , then 




Zq{{k + l)r) >1 — e-e “ 3 — 2/3r. 

a 

Proof. Recall the ODE for when Xc is not present is 

^ > ysVsiays - bys] - cys- 


By corollary 5.7 we know that for all t G [kr, {k + l)r], ys{t) > 2/o- The ODE for Zq is 

dzn 


dt 


klxaysZq - t 


Zq . 


(s,a,g)GA 


Since (s, a, q) G A, 


> k^XaysZq - k^XrZq 

> {ki - 5)xayoZq - {ks + 6)eZq 


During phase [kr + |, fcr + 3 ;^ > 1 _ so 


d z 

-^>{ki- (5)(1 - e)yoZq - (fcg + 6)eZq 
= aZq — /3Zq. 


By corollary 5.16 we have 


Zqikr + ^) > c{Zq) - - - e " 3 , 
6 a 
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and in the interval [kr + y, (A; + l)r], we have the inequalities 

dt - ' •> 

> -/3(l + 2e) 

> - 2 / 9 . 

Whence, at most 2/31 of Zg is lost during that time interval. Therefore, 

Zg{{k + l)r) > 1 — e — — — — 2/3r. 

a 


□ 


Lemma 5.24. Given construction 5.17 with X = Xa, let j/o and Zq be constants such 
that yo < E 2 from construction 5.21 If for every (s, a, g) G A, ysikr) < yo, and 
ZqikT) < Zq then 

Zq{{k + l)r) < Zo + 4d(fci + 5)yoT. 

Proof. Recall the ODE for yg when Xc is not present is 

^ < ysVsiays - bys] + cys- 


By corollary 5.11 we know that for all t G [kr, {k + l)r], ys{t) < yo- The ODE for Zq 

dzr. 


IS 


dt 


Therefore 


= ^ klXaysZq - klXrZq 

{s,a,q)eA 

A d[ki + (3)(1 + 2e)?/o(l + 2e) 
< 4:d{ki + 6)yo 


Zq{{k + l)r) < 2:0 + ^d{ki + 5)yQT 


□ 


Lemma 5.25. Given construction 5.17 with X = Xc, if Zq{kT) > zq, then 

where 


« = (^3 — ^)(1 — ^){zo — 2 e{k 3 + h)), and 

fd = (/c 3 + 5)(1 + 2e)(l + 2e + 2 e(/c 3 + 5) — ^ 0 ) T 4(A;4 + 6). 
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Proof. Recall the ODE for Zq is 


dzq 

dt 


kfXaVsZq - klXrZq 


(s,a,q)eA 


> —k^XrZ 


<? 


> —(^3 + “1“ 2e) 

> — 2e[k^ + 5 ). 


Therefore, at most 2e(/c3 + ^) y can be lost during the interval [kr, fcr + y 
for all t G [kr, /cr + y], 

Zg{t) > zo- 2e{kz + 5)t. 

The ODE for yq in the interval [kr + |, fcr + y] is 

^ = klXc[Zq% - ZqVq] + klVqVqlVq - Vq]] 

> (fcs - 5)(1 - e){zo - 2 e{k 3 + 6 ))yq 

— (fca + (5)(1 + 2e)(l + 2e + 2 e(/c 3 + 5) — zo)yq 

— (^4 + 5)(1 + 2 eYyq 

> ayq - /3yq. 


Now by lemma 5.14 we have 


2r. 


Vqi^T + y) > Vqikr + ^)e ^“+^^3 + c(Yq)(^ ^ ^ 


a 


)(1 




>c{Yq){ 


a 




(X P 

a 


)(i 


_ —(a+/3)J 


(X P 


)(i 




> 


a 


o T /S 


(1 


_ f: _ p“(“+4); 


)• 


Lemma 5.26. Given construction 5.17 with X = Xc, if Zqikr) < zq, then 

2 




(y-O+Oi y a) , 


where 


Q. — {k^ T (5)(1 T 2e)(^o T ^d{k\ T T 4(/c4 T 5), and 

(3 = {ks- 6){1 - e)(l - e - Zo - 4:d{ki + 6)e). 


. Thus, 


□ 
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Proof. Recall the ODE for Zq is 


dZq 

dt 


{s,a,q)&A 


klXaUsZq 


— k^x 


r^q 


< d{ki + (5)e(l + 2 e)(l + 2 e) 

< Adiki + 5)e 


Therefore, at most Ad{ki + 5)e‘^ can be gained dnring the interval [fcr, fcr + y]. Tims, 
for all t G [kr, fcr + y], 

Zq{t) < Zo + 4:d{ki + 5)eT. 

The ODE for i/q in the interval [kr + + y] is 

^ = klXc[Zqyq - ZqVq] + klVqVqlVq - Vq]] 

< {ks + 6){1 + 2e){zo + 4:d{ki + 6)eT)yq 

- {ks - S){1 - e){l - e - Zq - 4d{ki + 6)eT)yq 
+ (^4 + 5)(1 + 2e)‘^yq 

< ayq - /3yq. 


Now by lemma 5.14 we have 


2r, 


y) < ^)e + c(Fg)(^^^)(l - e 


a 


a + p 


a 


a 


= ciYq) ( + 

01 jj Oi jj 


3 —(q+/3) J 


= <y,) ( 




cr T /3 


^ — {a+fi)r 


3 -f 


a 


cr T /3 


< (1 + 2 e) 
2 




^-(a+/3) 


ycr T fd 


3 + 


a 


cr T /3 


< 


cr T /? 


(p- 


3 —(a+/3) J 


+ ofj 


□ 

Lemma 5.27. Let <h(ta) = ‘hi(w) and where $1 and $2 correspond to the 

first and second danse of the indnctive hypothesis <h. If constants e, S, y, r, 7 have 


25 












the constraints below for some e < 7 * < 7 , 


4(a + h) ’ 


1 . c < 


2 - "T > -r — log 


( 1 _ 7 *)(^ 3_1 + 7 ) 


+ T^log 


(1-7)(1-7*-E2) 


A\Ej, [(1- 7 )(E 3 - 1+ 7*)J E 2 [(1-7*)(1-7-E2)J/ ’ 

7> + 6 + e-(“+^)i, 

cr + p 


where the constants a, b, c, p, Ei, E 2 , E^, A are from construction 5.19 and where 
a = {ks- (5)(1 - e){zo - 2 e(fc 3 + 5)), 

/3 = (/c 3 + 5)(1 + 2e)(l + 2e + 2 e(/c 3 + ^) — ^ 0 ) “1“ 4(A;4 + <5), 

(^3 + b)e 


Zo = l-e- 


then 


{ki - 5 )(l-e)(l- 7 *) 


$(w) 


_ g-(fci-5)(i-p(i-T,*)r _ ^ 


$i(wa) 


Proof. Assume <h(tc) and let (s, a, q) G A. It suffices to show that |/q(3|tca|r) >1 — 7 . 
Let the constants a, h, c, p, Ei, E 2 , E 3 , A be as given in construction |5.19 , and let 
7 * be a constant such that 7 > 7 * > e. By lemma 5.20 and constraints (1) and (2), 

2/^((3|w;| + l)r) > 1 - 7 *. 


By lemma 5.23 and constraint (1) 

Zq{{3\w\ + 2)r) > zq, 


where zq is dehned above. By lemma 5.25 


//II \ 2 t , Q. 

9 ,(( 3 H+ 2 )r+-)> — 


e — e 


-(«+/3)5 


and by constraint (3) 


Finally, by corollary 5.7 


yqi{3\w\ + 2 )r+ y) > 1 - 7 . 


yq{3\wa\T) > 1 - 7 . 


□ 


Lemma 5.28. If the following constraints are satished with 7 > 7 * > El, 


1 / 1 . 

1 - r>— — log 


A* \El 


rb-Ei) 

7 ( 7 * - El) 


VW log 
E*2 


.T{El--f) 


2. 7 > 


Of “ 1 " 


{I3e 


-{a+yy. 


+ «)) 
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where 


01 — (/c3 + 5)(1 + 2e)(zo + 4:d{ki + 5))e + 4(^4 + <5), 

(3 = {ks - 6){1 - e)(l - e - zq - 4rf(fci + 5)e), 

^0 = 2 f + 2ci(fci + 6 )eT, 

\{k 3 - 6 ){l-e) ) 


then 


<h(w) $2(wa). 

Proof. Assume the constraints and <h(ta). Let q ^ 6 {I,wa). By lemma 


5.18 


^,((3|m| + 1 )t) < 2 ( + 2d(ki + 5) 

\{k 3 - 6 ){l-e) ) 


er 


= Zq. 


By lemma 5.22 and constraint (1), 

yq{i3\w\ + l)r) < 7 *. 


By lemma 5.24 we have 


By lemma 5.26 


2:q((3|M;| + 2)r) < Zq + Ad{ki + S)'y*T. 


yq{{3\w\ + 2)r + y) < + «) 

< 7- 


Finally, by corollary 5.11 


yq{3\wa\T) < 7 . 


□ 
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